#include <iostream>
#include <limits>
#include <cmath>
#include <algorithm>
using namespace std;

double eps = numeric_limits<double>::epsilon();

class Function
{
public:
    virtual double operator()(double _x) = 0;
    virtual double diff(double _x) {
        return 0;
    }
};
/*
class Qubic : Function
{
public:
    double operator()(double _x);
    double diff(double _x);
};
*/
class EquationSolver
{
public:
    virtual double solve() = 0;
};

//EquationSolver eq0;

//the bisection method.
class Bisection : public EquationSolver
{
private:
    double a;
    double b;
    double delta;
    int M;
    Function &f;
public:
    Bisection(double _a, double _b, double _delta, int _M, Function &_f) : a(_a), b(_b), delta(_delta), M(_M), f(_f){};
    double solve();
};
double Bisection::solve()
{
    double h = b - a;
    double u = f(a);
    double c;
    double w;
    for (int k = 1; k <= M; k++){
        h /= 2;
        c = a + h;
        w = f(c);
        if (fabs(h) < delta || fabs(w) < eps)
            break;
        else if (w * u >= 0)
            a = c;
    }
    return c;
}

//Newton's method.
class Newton : public EquationSolver
{
private:
    double x0, delta;
    int M;
    Function &f;
public:
    Newton(double _x0, double _delta, int _M, Function &_f) : x0(_x0), delta(_delta), M(_M), f(_f){};
    double solve();
};
double Newton::solve()
{
    double x = x0;
    double u;
    for (int k = 0; k <= M; k++){
        u = f(x);
        if (fabs(u) < eps)
            break;
        x = x - u / f.diff(x);
    }
    return x;
}

//the secant method.
class Secant : public EquationSolver
{
private:
    double x0, x1, delta;
    int M;
    Function &f;
public:
    Secant(double _x0, double _x1, double _delta, int _M, Function &_f) : x0(_x0), x1(_x1), delta(_delta), M(_M), f(_f){};
    double solve();
};

double Secant::solve()
{
    double xn = x1;
    double xn_1 = x0;
    double u = f(xn);
    double v = f(xn_1);
    double s;
    for (int k = 2; k <= M; k++)
    {
        if (fabs(u) > fabs(v)){
            swap(xn, xn_1);
            swap(u, v);
        }
        s = (xn - xn_1) / (u - v);
        xn_1 = xn;
        v = u;
        xn = xn - u * s;
        u = f(xn);
        if (fabs(xn - xn_1) < delta || fabs(u) < eps)
            break;
    }
    return xn;
}
